Supplement A: Estimating Travel Cost

Published

July 12, 2023

1 Overview

Goal: estimate cost-distance (in time) from each ranch to Ciudad Constitucion, La Paz, and local springs.

Data: a Digital Elevation Model (DEM) along with point locations for ranche clusters, springs, cities.

Method: Campbell’s hiking function (Campbell et al. 2019).

There are two watersheds in the project area. Each one has a single road that connects it to Highway 1 and thus to the markets in Ciudad Constitucion to the north and the capitol La Paz to the south. The “road” that connects the watersheds to each other is not maintained and thus rarely used. As a consequence, nearly all of the variance in travel time to cities is found within each watershed. For that reason, we estimate the time it takes to leave each watershed (“leaving” here means hitting the edge of the project area), rather than the true time it takes to get to each city. We interpret this as the relative difference in time it takes to get to Ciudad Constitucion from the northern watershed and the relative difference in time it takes to get to La Paz from the southern watershed. To get to La Paz from the northern watershed, we simply add an arbitrary half hour to the total time to get out of that watershed. The same would presumably be the case to get to Ciudad Constitucion from the southern watershed, though that value does not factor into our analysis.

For springs, Dr. MacFarlan’s ethnographic data suggest that individuals residing within the same ranch cluster tend to rely on one or two springs at most. Thus, we estimate road distance to the two nearest springs and take the average.

All spatial data are projected to the Mexico ITRF92 / UTM zone 12N CRS (EPSG:4485).

2 R Preamble

library(dplyr)
library(ggnewscale)
library(ggplot2)
library(here)
library(hiker)
library(sf)
library(terra)
library(viridis)

Specify consistent plot theme here:

Code
# this is designed to position the legend 
# at the top left corner of the plot area

theme_set(theme_bw(12))

theme_update(
  axis.text = element_text(size = rel(0.5), color = "gray45"),
  axis.text.y = element_text(angle = 90, hjust = 0.5),
  axis.ticks = element_line(color = "gray45", linewidth = 0.1),
  axis.ticks.length = grid::unit(0.1, "cm"),
  axis.title = element_blank(),
  legend.background = element_rect(fill = "transparent", color = "transparent"),
  legend.key.height = grid::unit(0.4, "cm"),
  legend.key.width = grid::unit(0.3, "cm"),
  legend.margin = margin(l = 0),
  legend.justification = c("left", "top"),
  legend.spacing.x = grid::unit(0.05, "cm"),
  legend.spacing.y = grid::unit(0.05, "cm"),
  legend.text = element_text(vjust = 0.7, margin = margin()),
  legend.title = element_blank(),
  panel.grid = element_blank(),
  plot.title = element_text(size = rel(1.1), face = "bold")
)

point_colors <- c("#FDE74C", "#EE7733", "#0077BB")

pretty_breaks <- function(x, .axis = "x", dx = 2500, n = 4){ 
  
  mn <- paste0(.axis, "min")
  mx <- paste0(.axis, "max")
  
  brks <- seq(x[[mn]] + dx, x[[mx]] - dx, length.out = n)
  
  round(brks, digits = -3)
  
}

3 Data

3.1 Geopackage

Specify path to geopackage database holding all spatial vector data.

gpkg <- here("data", "choyero.gpkg")

3.2 Project Window

bcs <- read_sf(gpkg, layer = "baja")

window <- read_sf(gpkg, layer = "window")

roads <- read_sf(gpkg, layer = "roads")

watersheds <- read_sf(gpkg, layer = "watersheds")

3.3 Elevation (DEM)

dem <- rast(here("data", "dem_30m.tif"))

3.4 O-D Points

# ORIGIN ---
clusters <- 
  read_sf(gpkg, layer = "clusters") |>
  mutate(variable = "cluster")

# DESTINATION --- to city, terminus at edge of window
terminus <- 
  read_sf(gpkg, layer = "terminus") |>
  mutate(variable = "terminus")

springs <- 
  read_sf(gpkg, layer = "springs") |>
  select(id, name) |>
  mutate(variable = "spring")
Code
od_points <- clusters |> 
  rename("name" = cluster) |> 
  mutate(name = as.character(name)) |> 
  select(name, variable) |> 
  bind_rows(terminus, springs) |> 
  mutate(variable = ifelse(variable == "terminus", "city", variable))

bb8 <- st_bbox(window)
dx <- bb8[["xmax"]] - bb8[["xmin"]]
dy <- bb8[["ymax"]] - bb8[["ymin"]]

# add hillshade relief to visualize elevation
# https://dominicroye.github.io/en/2022/hillshade-effects/
hillshade <- shade(
  slope = terrain(dem, "slope", unit = "radians"),
  aspect = terrain(dem, "aspect", unit = "radians"),
  angle = 45,
  direction = seq(45, 360, by = 45),
  normalize = TRUE
)

hillshade <- mask(sum(hillshade), vect(st_union(watersheds)))
names(hillshade) <- "shade"

gg <- ggplot() +
  geom_raster(
    data = as.data.frame(hillshade, xy = TRUE),
    aes(x, y, fill = shade)
  ) +
  scale_fill_distiller(
    palette = "Greys", 
    na.value = "transparent",
    guide = "none"
  ) +
  new_scale_fill() +
  geom_raster(
    data = as.data.frame(dem, xy = TRUE),
    aes(x, y, fill = elevation),
    alpha = 0.7
  ) +
  scale_fill_distiller(
    palette = "Greys", 
    na.value = "transparent",
    guide="none"
  ) +
  new_scale_fill() +
  geom_sf(
    data = watersheds, 
    fill = "transparent", 
    color = alpha("darkblue", 0.35),
    linewidth = 0.2
  ) +
  geom_sf(
    data = roads, 
    color = "black",
    linewidth = 0.3
  ) +
  geom_sf(
    data = od_points, 
    aes(shape = variable, fill = variable, color = variable),
    stroke = 0.3,
    size = 1.8
  ) +
  scale_color_manual(values = c("grey15", "white", "white")) +
  scale_fill_manual(values = alpha(point_colors, 0.75)) +
  scale_shape_manual(values = c(23, 22, 21)) +
  coord_sf(datum = 4485) +
  scale_x_continuous(
    breaks = pretty_breaks(bb8, "x"), 
    expand = expansion(add = 1000)
  ) +
  scale_y_continuous(
    breaks = pretty_breaks(bb8, "y", n = 3), 
    expand = expansion(add = 1000)
  ) +
  theme(legend.position = c(0.03, 0.98))

ggsave(
  plot = gg,
  filename = here("figures", "od-points.png"),
  dpi = 600,
  width = 3.5,
  height = 3.5 * (dy/dx)
)

4 Path analysis

To estimate least-cost paths, we first derive slope estimates from the DEM then calculate the hiking speed at which individuals can traverse any given cell using Campbell’s hiking function. After that, we decrease the cost in units of time to traverse grid cells crossed by roads and other paths using weights that approximate travel times by car or horse. This is equivalent to increasing the speed one travels in those cells. Below are the specific weights we use.

Down-weight for travel-cost
Type Weight Approx. Speed
primary 0.2 25 mph (41.5 kmh)
secondary 0.5 10 mph (16.2 kmh)
tertiary 0.8 07 mph (11.0 kmh)
primary_roads   <- roads |> filter(level == "P") |> st_buffer(dist = 30)
secondary_roads <- roads |> filter(level == "S") |> st_buffer(dist = 30)
tertiary_roads  <- roads |> filter(level == "T") |> st_buffer(dist = 30)

terrain <- dem |> 
  hf_terrain(max_slope = 35, decile = 30) |> 
  hf_channel(primary_roads,   .m = 0.2) |> 
  hf_channel(secondary_roads, .m = 0.5) |>
  hf_channel(tertiary_roads,  .m = 0.8)

remove(primary_roads, secondary_roads, tertiary_roads)
Code
gg <- ggplot() +
  geom_raster(
    data = as.data.frame(hf_rasterize(terrain), xy = TRUE),
    aes(x, y, fill = lyr.1)
  ) +
  scale_fill_viridis(
    name = "Seconds",
    breaks = c(10, 30, 50)
  ) +
  coord_sf(datum = 4485) +
  scale_x_continuous(
    breaks = pretty_breaks(bb8, "x"), 
    expand = expansion(add = 1000)
  ) +
  scale_y_continuous(
    breaks = pretty_breaks(bb8, "y"), 
    expand = expansion(add = 1000)
  ) +
  guides(
    fill = guide_colorbar(
      barwidth = 3.8, 
      barheight = 0.6,
      raster = FALSE,
      frame.linewidth = 0.2,
      ticks.linewidth = 0.7,
      title.position = "top"
    )
  ) +
  theme(
    legend.position = c(0.03, 0.98),
    legend.direction = "horizontal",
    legend.spacing.x = grid::unit(0.05, "cm"),
    legend.text = element_text(size = rel(0.5), margin = margin(l = 2)),
    legend.title = element_text(
      size = rel(0.8), 
      margin = margin(b = 2, t = 2),
      hjust = 0
    )
  )

ggsave(
  plot = gg,
  filename = here("figures", "cost-surface.png"),
  dpi = 600,
  width = 3.5,
  height = 3.5 * (dy/dx)
)

4.1 To springs

to_springs <- hf_hike(
  terrain, 
  from = clusters, 
  to = springs,
  add_cost = TRUE
)

to_springs <-
  to_springs |>
  group_by(from) |> 
  slice_min(cost, n = 2) |> # for each cluster, get two lowest cost paths
  mutate(cost = (cost/3600)) |> # convert to hours 
  summarize(
    cluster = unique(from),
    spring = paste(to, collapse = ","),
    cost = mean(cost)
  ) |> 
  select(cluster, spring, cost)

4.2 To cities

Please note that travel time is not strictly to cities, but to the edge of the project area. As noted above, for the northern watershed, the time to La Paz equals the time it takes to leave the watershed plus half an hour.

north <- clusters |> 
  st_filter(
     st_union(filter(watersheds, id %in% c(3294, 3296, 3299)))
  )

south <- clusters |> 
  st_filter(
    filter(watersheds, id == 3326)
  )

### NORTH ###
out_north <- hf_hike(
    terrain, 
    from = north, 
    to = terminus[1, ], 
    add_cost = TRUE
  ) |> 
  mutate(
    cluster = north$cluster,
    city = "Ciudad Constitucion", 
    cost = (cost/3600),
    watershed = "north"
  )

### SOUTH ###
out_south <- hf_hike(
    terrain, 
    from = south, 
    to = terminus[2, ], 
    add_cost = TRUE
  ) |> 
  mutate(
    cluster = south$cluster,
    city = "La Paz",
    cost = (cost/3600),
    watershed = "south"
  )

to_cities <- 
  bind_rows(out_north, out_south) |> 
  mutate(time_to_paz = ifelse(watershed == "north", cost + 0.5, cost)) |>
  select(cluster, city, cost, time_to_paz)

remove(north, south, out_north, out_south)

4.3 LCP Map

Code
paths <- bind_rows(
  to_springs |> mutate(destination = "To springs"),
  to_cities |> mutate(destination = "To cities")
)

gg <- ggplot() +
  geom_raster(
    data = as.data.frame(hillshade, xy = TRUE),
    aes(x, y, fill = shade)
  ) +
  scale_fill_distiller(
    palette = "Greys", 
    na.value = "transparent",
    guide = "none"
  ) +
  new_scale_fill() +
  geom_raster(
    data = as.data.frame(dem, xy = TRUE),
    aes(x, y, fill = elevation),
    alpha = 0.7
  ) +
  scale_fill_distiller(
    palette = "Greys", 
    na.value = "transparent",
    guide="none"
  ) +
  new_scale_fill() +
  geom_sf(
    data = watersheds, 
    fill = "transparent", 
    color = alpha("darkblue", 0.35),
    linewidth = 0.2
  ) +
  geom_sf(
    data = roads, 
    color = "black",
    linewidth = 0.3
  ) +
  geom_sf(
    data = paths |> st_buffer(125) |> group_by(destination) |> summarize(),
    fill = "#C7FFED",
    color = "#00291D",
    linewidth = 0.2
  ) +
  facet_wrap(vars(destination)) +
  geom_sf(
    data = od_points, 
    aes(shape = variable, fill = variable, color = variable),
    stroke = 0.3,
    size = ifelse(od_points$variable == "city", 2.2, 1.5) |> rep(2)
  ) +
  scale_color_manual(values = c("grey15", "white", "white")) +
  scale_fill_manual(values = alpha(c(point_colors, 0.75))) +
  scale_shape_manual(values = c(23, 22, 21)) +
  coord_sf(datum = 4485) +
  scale_x_continuous(
    breaks = pretty_breaks(bb8, "x"), 
    expand = expansion(add = 1000)
  ) +
  scale_y_continuous(
    breaks = pretty_breaks(bb8, "y", n = 3), 
    expand = expansion(add = 1000)
  ) +
  theme(
    legend.position = c(0.01, 0.98),
    strip.background = element_blank(),
    strip.text = element_text(size = rel(1.1), hjust = 0)
  ) 

ggsave(
  plot = gg,
  filename = here("figures", "least-cost-paths.png"),
  dpi = 600, 
  width = 5.75,
  height = 5.75 * (dy/(2*dx)) + 0.47
)

4.4 Save

write_sf(
  to_springs,
  dsn = gpkg,
  layer = "path_to_springs"
)

write_sf(
  to_cities,
  dsn = gpkg,
  layer = "path_to_cities"
)

5 References

Campbell, Michael J., Philip E. Dennison, Bret W. Butler, and Wesley G. Page. 2019. “Using Crowdsourced Fitness Tracker Data to Model the Relationship Between Slope and Travel Rates.” Applied Geography 106: 93–107. https://doi.org/https://doi.org/10.1016/j.apgeog.2019.03.008.

6 Session Info

Code
# save the session info as an object
pkg_sesh <- sessioninfo::session_info(pkgs = "attached")

# inject the quarto info
pkg_sesh$platform$quarto <- paste(
  quarto::quarto_version(), 
  "@", 
  quarto::quarto_path()
  )

# print it out
pkg_sesh
─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.3.1 (2023-06-16 ucrt)
 os       Windows 11 x64 (build 22621)
 system   x86_64, mingw32
 ui       RTerm
 language (EN)
 collate  English_United States.utf8
 ctype    English_United States.utf8
 tz       America/Denver
 date     2023-07-12
 pandoc   3.1.1 @ C:/Program Files/RStudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)
 quarto   1.4.176 @ C:\\PROGRA~1\\Quarto\\bin\\quarto.exe

─ Packages ───────────────────────────────────────────────────────────────────
 package     * version date (UTC) lib source
 dplyr       * 1.1.2   2023-04-20 [1] CRAN (R 4.3.1)
 ggnewscale  * 0.4.9   2023-05-25 [1] CRAN (R 4.3.1)
 ggplot2     * 3.4.2   2023-04-03 [1] CRAN (R 4.3.1)
 here        * 1.0.1   2020-12-13 [1] CRAN (R 4.3.1)
 hiker       * 0.1.0   2023-07-12 [1] local
 sf          * 1.0-13  2023-05-24 [1] CRAN (R 4.3.1)
 terra       * 1.7-39  2023-06-23 [1] CRAN (R 4.3.1)
 viridis     * 0.6.3   2023-05-03 [1] CRAN (R 4.3.1)
 viridisLite * 0.4.2   2023-05-02 [1] CRAN (R 4.3.1)

 [1] C:/Users/kenne/AppData/Local/R/win-library/4.3
 [2] C:/Program Files/R/R-4.3.1/library

──────────────────────────────────────────────────────────────────────────────